Linear scaling calculation of a n-type GaAs quantum dot 
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A linear scale method for calculating electronic properties of large and complex systems is intro- 
duced within a local density approximation. The method is based on the Chebyshev polynomial 
expansion and the time-dependent method, which is tested in calculating the electronic structure of 
a model n-type GaAs quantum dot. 
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I. INTRODUCTION 



Linear scale methods for calculating the electronic 
structures have been actively investigated in the last 
decade because of increasing demands for predicting 
properties of large and complex systems with computa- 
tional cost linear scale with respect to the system size 
N [l|. There are several approaches for achieving linear 
scaling, such as the divide-and-conquer (DC) method 3j 
the density-matrix minimization (DMM) method Q , the 
orbital minimization (OM) method 0], and the Cheby- 
shev polynomial expansion (CPE) method Com- 
putational efficiency and applicability for specific systems 
have been mostly tested based on the tight-binding (TB) 
formalism. The DC method divides a system into subsys- 
tems in physical space and obtains the density matrix for 
each subsystem. This method is highly efficient if small 
localization region can be chosen as subsystems, but this 
depends on the problems and becomes more difficult for 
calculations based on the finite-difference (FD) formal- 
ism with a large basis set. While the TB method is very 
successful in quantum chemistry, care must be taken for 
constructing appropriate basis set for a particular prob- 
lem 0. A calculation based on the FD formalism Q is 
straightforward and is widely used for electronic struc- 
ture calculations of semiconductors and biochemical sys- 
tems. The DMM and OM methods, which require to 
store the whole density matrix and the whole Wannier 
functions, respectively, suffer from their large memory 
requirements. In the CPE method, the memory require- 
ments are significantly reduced because only small num- 
ber of column vectors is required to be stored. Since 
neither division into subsystems or the initial guess of 
the initial state is required, the CPE method is straight- 
forwardly applied to a wide variety of systems. The 
other important advantage of the CPE method is suit- 
ability for parallel implementation. Because the most 
time-consuming part of the calculation is matrix-times- 
vector multiplication, where each column of the Hamilto- 
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nian matrix can be treated as independent, communica- 
tions between clusters are minimized. The CPE method 
is thus suitable for achieving linear scaling based on the 
FD formalism with a large basis set. 

In the CPE method the electron density is evaluated 
by using a matrix representation of the Fermi-operator, 
which is expanded in the Chebyshev matrix polynomi- 
als. The so-called Gibbs oscillation in the zero tem- 
perature case is suppressed by using finite-temperature 
Fermi operator. |T3| In the tight binding approach, 
the linear scaling is obtained by a truncated Hamilto- 
nian which retains only matrix elements inside a local- 
ization region Reasonably small localization can 
be defined for a tight-binding approach with, for exam- 
ple, atom-centered basis functions. In the FD formal- 
ism, it is not obvious how to define a localization region 
where basis functions are retained. Moreover, because 
the number of basis functions within a localization region 
becomes much larger than the tight binding approach, 
the crossover point where the linear scaling approach is 
faster than a conventional approach such as a conjugate 
gradient method (CGM) becomes significantly larger. 

This leads us to utilize the other approach of calculat- 
ing the trace of a large matrix by using random vectors 
calculating physical quantities such as energy, 
electron density, or linear response function, the trace of 
a relevant operator A needs to be calculated. If A is ex- 
pressed in terms of a basis set (j)q,q = 1, Nd as tT[A] = 
J2q=i ^^qq' calculation of this part costs 0{N^) if the 
matrix is expanded in the Chebyshev matrix polynomi- 
als [l| . By introducing a random phase vector as defined 

by 1$) = J2q=i \l)^q^ where {\q)} is a basis set and £,q 
are a set of random phase variables, the trace is evalu- 
ated at the cost of 0{Nd) as given by tr[^] = {{$\A\$)), 
where ((•)) stands for statistical average. The overall 
linear scaling is obtained by this method. The random 
phase vector was shown to give results with the smallest 
statistical error [l^. This approach is also known to 
show a useful feature called the self-averaging effect that 
the ffuctuation in some physical quantities decreases with 
increase in Nd for sparse or banded matrices Anm ■ With 
a combination of this approach and the time-dependent 
method di] (CPE-TDM), Hnear response functions or 
electron density of states (DOS) are calculated by inte- 
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grating the time-dependent Schrodinger equation with- 
out calculating eigenenergies or eigenstates. The compu- 
tational time of CPE-TDM scales as 0{N), as compared 
with that of the conventional method such as conjugate 
gradient method (CGM), which grows as 0{N^). Thus 
CPE-TDM enables us to calculate electronic properties of 
large systems which require prohibitively large computa- 
tional time by CGM. CPE-TDM was applied to calculate 
the optical properties of hydrogenated Si nanocrystals 
containing atoms more than 10,000 within the empirical 
pseudopotential formalism [l^, [3l , the optical properties 
of carbon nanocrystals and polysilane [ll] , and the 
electron spin resonance spectrum of s = 1/2 antiferro- 
magnet Cu benzoate jl^ , which have proved the advan- 
tages of CPE-TDM. However, CPE-TDM has not been 
applied to calculation of the electronic structure within 
a local density approximation (LDA). Applications of a 
linear scaling method with the self-consistent-field level 
of theory are still very limited, but this level of calcula- 
tion using Gaussian basis sets has been demonstrated to 
be practical, [l^l 

In this paper, we report on an implementation of 
CPE-TDM for a large scale calculation of the electronic 
structure of n-type GaAs quantum dot (QD) [U [S^l 
within a LDA based on a FD formalism and compare 
the results with a CGM. 



II. METHOD 

The model structure is a 20 nm-wide GaAs quantum 
well sandwiched by undoped Alj;Gai_3;As {x — 0.3) bar- 
riers, which confine the electrons with the effective-mass 
TO* in the z direction. For QDs, the electrons are assumed 
to be laterally confined to a harmonic oscillator with fre- 
quency cjq, which may be created by a surface gate struc- 
ture j22| in experiments. The electrons are assumed to 
be supplied from 5 nm-thick Si-doped AlxGai_xAs layer, 
located 20 nm above the GaAs quantum well layer. The 
Fermi-energy (Ep) is taken as the origin of the energy. 
The Fermi-level pinning model is assumed [23| . The num- 
ber of the electrons in a QD is not fixed to an integer 
number and is determined by Ep and the potential en- 
ergy. 

The model Hamiltonian of the system within the LDA 



H = ^ + lrn*cuUx^ + 2/2) + V, (z) + Vh (r) + (r) 

(1) 

where Vc{z), Vf/(r), and Vx{r) are the vertical confining 
potential, the Hartree potential, and the exchange poten- 
tial, respectively. A 3D mesh of 64 x 64 x 8 is used for the 
calculation of the electron density, and 64 x 64 x 16 is used 
for the calculation of the potentials. The axis perpendic- 
ular to the quantum well layer is taken to be z-direction 
and the grid-spacing Ax is fixed to be 5 nm. The Hamil- 



tonian is discretized in real-space by the higher-order fi- 
nite difference method (2^ . [25| . 

The electron density at finite temperature is given by 



n{v) ^Y.^*{v)cly,{v)}{{E, - Ep)) 



(2) 



where (f)j and Ej are the one-particle wave function and 
the energy of the jth electron state, respectively, which 
are obtained by CGM. f{E.j — Ep) ~ ^Ej-Ep) ^ is the 
Fermi distribution function at inverse temperature (3. We 
use (3 = 4000 eV~^ corresponding to the temperature 
T — 2.9 K. The electron states above Ep are partially 
occupied due to this finite temperature effect. The in- 
troduction of finite temperature accelerate convergence 
of the self-consistent-field loop. 

In CPE-TDM, a random phase vector as defined by 
1^) = \l)iq^ where are a set of random phase 



variables 



is used as an initial state. Here $ is 



a Nx X Ny X Nz column vector for a system defined by a 
real-space uniform grid of x Ny x N^- The electron 
density n{v) is extracted by the Fermi operator function 
f{H) 



g/3(H-i3p)_|_]^ 



as 



n(r) = ((|($|/(i7)|r)p)). 



(3) 



where j3 is connected to a real temperature. The Fermi 
operator is evaluated by the Chebyshev polynomial ex- 
pansion, 



/(i/)|$>=^afc(/3)n(F)|$). 



(4) 



The length of the Chebyshev expansion for precision 
is given 



10--° is given by [i,[T3] 



P 



-iD-l)(3-AE 



(5) 



where AE = {E^ax~ E^,n) /2. We use D ^ &, AE ^ 1.0 
eV, 13 = 4000 eV-\ giving P = 13333. A calculation 
was also performed with I? = 9, and we find that the 
differences in the total number of electrons (A^e) and the 
Hartree potential {Vh) between the two cases of _D = 6 
and Z) = 9 were less than 1 x 10"'^ and 1 x 10^^ eV, 
respectively. 

The electron density is calculated with jmax sets of |<I>) 
as n(r) = E^Ii"'"' {AfWWj) /jmax- The fiuctuation for 
the random phase vector is [H, H^l 



5H/L ■ 



V2 



2m* [Ax)^ V>J' 



(6) 



where L = Nx{Ny)h as the number of meshes N ^ oo. 
The statistical error decreases as 1 / \Jjmax in general. 

While it is known that other representation of a 
smoothed step function such as a complementary error 
function yields improvements of degree of polynomial ex- 
pansion [28| . we use the Fermi operator because this 



is physically correct for electronic structure calculations 
at finite temperature. The Hartree and exchange po- 
tentials are calculated by using Eq. (3). Therefore, it 
is not necessary to obtain eigenvalues or eigenfunctions. 
The new solution of the potential V^°^{r) is combined 
with the solution obtained for the previous iteration by 
VH(r) = (1 - a)V'^'^{r) + a Vg<'^(r). Similarly, in order 
to reduce the statistical fluctuation, n{r) is combined 
with the density obtained for the previous iteration by 
n(r) = (1 - 7)n°''^(r) + 7n"™(r). The parameter a is 
fixed to be 0.08, and the parameter 7 is varied between 
0.3 and O.f. 

A real-time Green's function G{lu£ + irj) is calculated 
by a time evolution method by solving a homogeneous 
Schrodinger equation numerically with an initial condi- 
tion (j){q,t = 0) = \q) as ^ 

Mq> T) = i-i) r dt'c^iq, t')e'('-^+'^)^' (7) 
Jo 



^G{LUi+lTj)\q). (9) 

This method is as efficient as the CPE method with a 
carefully chosen Gibbs damping factor. [6, 9] The DOS 
is then calculated at the cost of 0{Nd) as given by 




p{Lu) = ~-lin{TT[G{Lu + ir])]), (10) 



= -llm(^((e'(*'-V)))(g|G(c. + t^)\q')) (11) 

9,9 

= --Im((( ($|G(c. + *,7)|$) ))). (12) 

TT 

The DOS is calculated with fcmax sets of |<i>). The 
energy resolution 77 is chosen to be 0.25 meV. It should 
be noted that /cmax used for calculating the DOS can be 
independently chosen from jmax for each self-consistent 
iteration procedures. 



III. RESULTS AND DISCUSSIONS 

Model calculations are performed for GaAs QDs con- 
taining about 77 electrons. We take uiq — 3 meV 
for a typical GaAs QD [2lj. The number of the self- 
consistent iterations is fixed to 100 for both the CGM 
and CPE-TDM calculations. The potential is converged 
to \Vii{r) - V^'='"{r)\ < 0.003 meV for the CGM cal- 
culation. The electron density distributions are shown 



FIG. 1: (a) The electron density distribution obtained by 
CPE-TDM. 128 sets of random vectors are used at each self- 
consistent iteration procedure, (b) The electron density dis- 
tribution obtained by CGM. 



in Fig. [U for CPE-TDM with j^ax = 128 and CGM. 
The calculated electron density distribution reasonably 
agrees with the result by a CGM within the statistical 
fluctuations. The Friedel-type spatial oscillations of the 
electron density [2^ are reproduced in both the results 
by the CPE-TDM and CGM. 

The calculated Hartree potentials reasonably agree 
with the potential obtained by the CGM as shown in 
Fig. [5] (a). Differences of the calculated Hartree poten- 
tials with that by the CGM are examined in Fig. [2] (b). 
The absolute values of the difference are smaller than 1.0 
and 0.4 meV for jmax = 8 and 16, respectively. Figure [2] 
(c) shows that the standard deviations of the differences 
of the calculated Hartree potentials follows the curve pro- 
portional to 1/ -y/jinax as expected. 

The calculated DOS are shown in Fig. E] For CPE- 
TDM, the self-consistent iteration procedures are per- 
formed with jmax — 8, 16, 32, 64, and 128. The same 
number of random phase vectors are used for evaluating 
p{llj) except for the case of jmax = 8 where p{uj) is eval- 
uated with fcinax = 8 and 64. It can be seen that the 
statistical fluctuations decrease with increase in jmax in 
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FIG. 2: (a) Cross-sectional views of the calculated Hartree 
potentials on the plane at the center of the quantum well 
layer obtained by CPE- TDM (^^^'^-'^^^(r)) with (i) 8, (ii) 
16, (iii) 32, (iv) 64, and (v) 128 sets of random phase vectors 
for extracting n(r) at each self consistent iteration procedure, 
and (vi) Vh{^) obtained by CGM. (b) Differences of obtained 
Hartree potentials V^^'^-'^°'^ [r) - V§^^{r) with (i) 8, (ii) 
16, (iii) 32, (iv) 64, and (v) 128 sets of random phase vectors, 
(c) Standard deviations of the calculated Hartree potentials 
depending on the number of random phase vectors at each 
self consistent iteration procedure jmax. The best fitted curve 
proportional to l/\/jmax is also shown. 



calculating pitu). There are two types of the fluctuations 
observed in Fig. [S] One is the fluctuation in the peak 
energy positions, and the other is the fluctuation in the 
peak heights. The former can be reduced by increasing 
Jmax and by decreasing the mixing parameter 7. The 
latter also depends on fcmax- In fact, the fluctuations in 
the peak heights are reduced by increasing fcmax from 8 
to 64 with small changes in the peak energy positions in 
the case of jmax = 8 as shown in Fig. [3] (b). Figure [3] (b) 
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FIG. 3: (a) Density of states {p(uj)) obtained by CPE-TDM 
with jmax = fcnax = (i) 8, (ii) 16, (iii) 32, (iv) 64, and (v) 128 
sets of random phase vectors for extracting n(r) at each self 
consistent iteration procedure and for evaluating p{u!). (vi) 
p(a;) obtained by CPE-TDM with 8 sets of random phase vec- 
tors for extracting n(r) and fcmax = 64 sets of random phase 
vectors for evaluating p{uj). (vii) p(w) obtained by CGM. 
(b) Standard deviations of the difference of the peak heights 
of the DOS obtained by CPE- TDM and by CGM depend- 
ing on Jmax (solid circles). Standard deviation for CPE-TDM 
with jmax = 8 at each self consistent iteration procedure and 
fcmax=64 for evaluating p{uj) is shown (open circle) The best 
fitted curve proportional to 1 / V jmax is also shown, (c) Total 
number of electrons (A^e) depending on jmax- The horizontal 
line shows A^e — 77.1 obtained by CGM. 



shows that the standard deviations of the peak heights 
also follows the curve proportional to 1/ V fcmax . 

Finally we note that the statistical fluctuation of the 
total number of electrons (TVg) is smaller than that of 
DOS because of the self-averaging effect. Figure [3] (c) 
shows calculated N,, depending on jmax- The statisti- 
cal errors as compared with Ne = 77.1 by CGM are as 
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small as 2% for jmax = 8, which indicates that the self- 
averaging effect is effective for a sparse banded matrix 
case as illustrated in this paper. 

Our linear scale method opens up possibilities for cal- 
culating the electronic and optical properties of large and 
complex systems, such as QD arrays with interaction be- 
tween QDs and devices employing the Rashba type spin- 
orbit interaction [s^. It should also be possible to cal- 
culate the electronic structure of nanostructures within 
a LDA with ab initio pseudopotentials. Because the 



Green's function can be effectively estimated by CPE- 
TDM, the properties of the electronic system such as the 
DC and Hall conductivities, and the optical absorption 
spectra, are obtained within 0{N) computational costs. 

In conclusions, it has been demonstrated that CPE- 
TDM can be applied to a large scale calculation of a 
model QD within a LDA based on a FD formalism de- 
spite the presence of the statistical fluctuations of the 
calculated quantities originated from the random phase 
vectors. 
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